% Cell Cluster Analysis
% Eden Ford
% 30 May 2023

% Code to calculate (1) the number of clusters per volume and (2) the
% number of cells per cluster

% General definitions:
% Object - cell or cluster of cells that are touching/interacting with one
% another

close all
clear all
clc

%% 0 mM CMP1a Data Import
filename = 'Cells.xlsx'; % file name to read
sheet = 1; % sheet to read

[num_0,txt_0,raw_0] = xlsread(filename,sheet);

ImageNames_0 = raw_0(2:end,2); % list of image each object is found in
ImageNumber_0 = unique(ImageNames_0); % create array of unique image names
ImageNamesChar_0 = char(ImageNames_0); % convert cell aray into character array
ImageNumberChar_0 = char(ImageNumber_0); % convert cell aray into character array
Population_0 = raw_0(2:end,3); % list of population type
PopulationChar_0 = char(Population_0); % list of population type
Volume_0 = cell2mat(raw_0(2:end,7)); % list of shape factor of each object
ShapeFactor_0 = cell2mat(raw_0(2:end,6)); % list of shape factor of each object

T_Full_0 = table(ImageNames_0,Population_0,Volume_0,ShapeFactor_0) % organize all data into table

%% 10 mM CMP1a Data Import
filename = 'Cells.xlsx'; % file name to read
sheet = 2; % sheet to read

[num_10,txt_10,raw_10] = xlsread(filename,sheet);

ImageNames_10 = raw_10(2:end,2); % list of image each object is found in
ImageNumber_10 = unique(ImageNames_10); % create array of unique image names
ImageNamesChar_10 = char(ImageNames_10); % convert cell aray into character array
ImageNumberChar_10 = char(ImageNumber_10); % convert cell aray into character array
Population_10 = raw_10(2:end,3); % list of population type
PopulationChar_10 = char(Population_10); % list of population type
Volume_10 = cell2mat(raw_10(2:end,7)); % list of shape factor of each object
ShapeFactor_10 = cell2mat(raw_10(2:end,6)); % list of shape factor of each object

T_Full_10 = table(ImageNames_10,Population_10,Volume_10,ShapeFactor_10); % organize all data into table

%% Cell clusters per 3 images

% Initialize cell cluster counts
ClusterCount_0 = [0 0 0]; % cells cluster count in 0 mM [gel1 gel2 gel3]
ClusterCount_10 = [0 0 0]; % cells cluster count in 10 mM [gel1 gel2 gel3]

for i = 1:length(ImageNamesChar_0)
    if ImageNamesChar_0(i,1:11) == '0mMCMP_Gel1'
        if PopulationChar_0(i,1:4) == 'F-Ac'
            ClusterCount_0(1) = ClusterCount_0(1) +1;
        end
    elseif ImageNamesChar_0(i,1:11) == '0mMCMP_Gel2'
        if PopulationChar_0(i,1:4) == 'F-Ac'
            ClusterCount_0(2) = ClusterCount_0(2) +1;
        end
    elseif ImageNamesChar_0(i,1:11) == '0mMCMP_Gel3'
        if PopulationChar_0(i,1:4) == 'F-Ac'
            ClusterCount_0(3) = ClusterCount_0(3) +1;
        end
    end
end

for i = 1:length(ImageNamesChar_10)
    if ImageNamesChar_10(i,1:12) == '10mMCMP_Gel1'
        if PopulationChar_10(i,1:4) == 'F-Ac'
            ClusterCount_10(1) = ClusterCount_10(1) +1;
        end
    elseif ImageNamesChar_10(i,1:12) == '10mMCMP_Gel2'
        if PopulationChar_10(i,1:4) == 'F-Ac'
            ClusterCount_10(2) = ClusterCount_10(2) +1;
        end
    elseif ImageNamesChar_10(i,1:12) == '10mMCMP_Gel3'
        if PopulationChar_10(i,1:4) == 'F-Ac'
            ClusterCount_10(3) = ClusterCount_10(3) +1;
        end
    end
end
ClusterCount_0
ClusterCount_10

% ANOVA
ClusterANOVA = [ClusterCount_0',ClusterCount_10'];
[p_Cluster tbl_Cluster stats_Cluster] = anova1(ClusterANOVA) % where T is a matrix with data columns
Cluster_comparison = multcompare(stats_Cluster) % multiple comparisons test

%% Cells per cluster

% Initialize cell cluster counts
CellPerClusterCount_0_Gel1 = []; % cells per cluster count in 0 mM
CellPerClusterCount_0_Gel2 = []; % cells per cluster count in 0 mM
CellPerClusterCount_0_Gel3 = []; % cells per cluster count in 0 mM
CellPerClusterCount_10_Gel1 = []; % cells per cluster count in 10 mM
CellPerClusterCount_10_Gel2 = []; % cells per cluster count in 10 mM
CellPerClusterCount_10_Gel3 = []; % cells per cluster count in 10 mM

for i = 1:length(ImageNamesChar_0)
    if ImageNamesChar_0(i,1:11) == '0mMCMP_Gel1'
        if PopulationChar_0(i,1:4) == 'F-Ac'
            CellPerClusterCount_0_Gel1 = [CellPerClusterCount_0_Gel1,0];
        elseif PopulationChar_0(i,1:4) == 'DAPI'
            CellPerClusterCount_0_Gel1(end) = CellPerClusterCount_0_Gel1(end)+1;
        end
    elseif ImageNamesChar_0(i,1:11) == '0mMCMP_Gel2'
        if PopulationChar_0(i,1:4) == 'F-Ac'
            CellPerClusterCount_0_Gel2 = [CellPerClusterCount_0_Gel2,0];
        elseif PopulationChar_0(i,1:4) == 'DAPI'
            CellPerClusterCount_0_Gel2(end) = CellPerClusterCount_0_Gel2(end)+1;
        end
    elseif ImageNamesChar_0(i,1:11) == '0mMCMP_Gel3'
        if PopulationChar_0(i,1:4) == 'F-Ac'
            CellPerClusterCount_0_Gel3 = [CellPerClusterCount_0_Gel3,0];
        elseif PopulationChar_0(i,1:4) == 'DAPI'
            CellPerClusterCount_0_Gel3(end) = CellPerClusterCount_0_Gel3(end)+1;
        end
    end
end

for i = 1:length(CellPerClusterCount_0_Gel1)
    if CellPerClusterCount_0_Gel1(i) == 0
        CellPerClusterCount_0_Gel1(i) = 1;
    end
end

for i = 1:length(CellPerClusterCount_0_Gel2)
    if CellPerClusterCount_0_Gel2(i) == 0
        CellPerClusterCount_0_Gel2(i) = 1;
    end
end

for i = 1:length(CellPerClusterCount_0_Gel3)
    if CellPerClusterCount_0_Gel3(i) == 0
        CellPerClusterCount_0_Gel3(i) = 1;
    end
end
CellPerClusterCount_0 = [CellPerClusterCount_0_Gel1 CellPerClusterCount_0_Gel2 CellPerClusterCount_0_Gel3]';
length(CellPerClusterCount_0);

for i = 1:length(ImageNamesChar_10)
    if ImageNamesChar_10(i,1:12) == '10mMCMP_Gel1'
        if PopulationChar_10(i,1:4) == 'F-Ac'
            CellPerClusterCount_10_Gel1 = [CellPerClusterCount_10_Gel1,0];
        elseif PopulationChar_10(i,1:4) == 'DAPI'
            CellPerClusterCount_10_Gel1(end) = CellPerClusterCount_10_Gel1(end)+1;
        end
    elseif ImageNamesChar_10(i,1:12) == '10mMCMP_Gel2'
        if PopulationChar_10(i,1:4) == 'F-Ac'
            CellPerClusterCount_10_Gel2 = [CellPerClusterCount_10_Gel2,0];
        elseif PopulationChar_10(i,1:4) == 'DAPI'
            CellPerClusterCount_10_Gel2(end) = CellPerClusterCount_10_Gel2(end)+1;
        end
    elseif ImageNamesChar_10(i,1:12) == '10mMCMP_Gel3'
        if PopulationChar_10(i,1:4) == 'F-Ac'
            CellPerClusterCount_10_Gel3 = [CellPerClusterCount_10_Gel3,0];
        elseif PopulationChar_10(i,1:4) == 'DAPI'
            CellPerClusterCount_10_Gel3(end) = CellPerClusterCount_10_Gel3(end)+1;
        end
    end
end

for i = 1:length(CellPerClusterCount_10_Gel1)
    if CellPerClusterCount_10_Gel1(i) == 0
        CellPerClusterCount_10_Gel1(i) = 1;
    end
end

for i = 1:length(CellPerClusterCount_10_Gel2)
    if CellPerClusterCount_10_Gel2(i) == 0
        CellPerClusterCount_10_Gel2(i) = 1;
    end
end

for i = 1:length(CellPerClusterCount_10_Gel3)
    if CellPerClusterCount_10_Gel3(i) == 0
        CellPerClusterCount_10_Gel3(i) = 1;
    end
end
CellPerClusterCount_10 = [CellPerClusterCount_10_Gel1 CellPerClusterCount_10_Gel2 CellPerClusterCount_10_Gel3]';
length(CellPerClusterCount_10);

%% Cells per Cluster Distribution and Statistical Analysis

%ShapeFactor_0
%ShapeFactor_10

% Display results
disp(' ');
disp('cells per cluster analysis');
disp(' ');

% Moments of distribution
Type_ShapeFactor = {'0mMCMP'; '10mMCMP'};
% MEAN
Average_ShapeFactor = [mean(CellPerClusterCount_0); mean(CellPerClusterCount_10)];
% STANDARD ERROR (DERIVED FROM VARIANCE)
StError_ShapeFactor = [std(CellPerClusterCount_0/sqrt(3)); std(CellPerClusterCount_10/sqrt(3))];
% MEDIAN
Median_ShapeFactor = [median(CellPerClusterCount_0); median(CellPerClusterCount_10)];
% MODE
Mode_ShapeFactor = [mode(CellPerClusterCount_0); mode(CellPerClusterCount_10)];
% SKEWNESS
Skewness_ShapeFactor = [skewness(CellPerClusterCount_0); skewness(CellPerClusterCount_10)];
% KURTOSIS
Kurtosis_ShapeFactor = [kurtosis(CellPerClusterCount_0); kurtosis(CellPerClusterCount_10)];

ShapeFactorTable = table(Type_ShapeFactor, Average_ShapeFactor, StError_ShapeFactor, Median_ShapeFactor, Mode_ShapeFactor, Skewness_ShapeFactor, Kurtosis_ShapeFactor);
ShapeFactorTable.Properties.VariableNames = {'GelComposition' 'Average' 'StError' 'Median' 'Mode' 'Skewness' 'Kurtosis'};
disp(ShapeFactorTable);

% Statistical analysis
filename = 'Cells.xlsx'; % file name to read
sheet = 3; % sheet to read

[num,txt,raw] = xlsread(filename,sheet);
CellsPerClusterFull = num(3:end,:);

[p tbl stats] = anova1(CellsPerClusterFull) % where T is a matrix with data columns
CellsPerCluster_comparison = multcompare(stats) % multiple comparisons test[h2a_RL3,p2a_RL3] = ttest2(ShapeFactor_0mMCMP,ShapeFactor2a_Low);

%% Histograms

figure
BinWidth = 2;
xMax = 70;
yMax = 1;
edges = linspace(0, 1, 21);
histogram_distribution = zeros(6,20);

    % 0 mM CMP
    H2a = subplot(2,1,1);
    hold on
    histogram(CellPerClusterCount_0, 'Normalization', 'probability', 'BinWidth', BinWidth);
    ylabel('Fraction of Objects')
    hold off
    title('0 mM CMP')
    xlim([0 xMax])
    ylim([0 yMax])
    histogram_distribution(1,:) = histcounts(CellPerClusterCount_0,edges);
    
    % 10 mM CMP
    H2b = subplot(2,1,2);
    hold on
    histogram(CellPerClusterCount_10, 'Normalization', 'probability', 'BinWidth', BinWidth);
    % ylabel('Fraction of Objects')
    hold off
    title('10 mM CMP')
    xlim([0 xMax])
    ylim([0 yMax])
    histogram_distribution(2,:) = histcounts(CellPerClusterCount_10,edges);

